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Abstract 

The reweighting method is widely used in numerical studies of QCD, in particular, for the cases in 
which the conventional Monte-Carlo method cannot be applied directly, e.g., finite density QCD. However, 
the application range of the reweighing method is restricted due to several problems. One of the most 
severe problems here is the overlap problem. To solve it, we examine a multipoint reweighting method 
in which simulations at several simulation points are combined in the data analyses. We systematically 
study the applicability and limitation of the multipoint reweighting method in two-flavor QCD at zero 
density. Measuring histograms of physical quantities at a series of simulation points, we apply the multipoint 
reweighting method to calculate the meson masses as continuous functions of the gauge coupling /3 and the 
hopping parameters k. We then determine lines of constant physics and beta functions, which are needed 
in a calculation of the equation of state at finite temperature. 
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I. INTRODUCTION 


At extremely high temperatures and/or densities, the quark matter is expected to turn into new 
phases. Clarification of the nature of these states as well as the phase structure of QCD is important 
in understanding the evolution of the Universe around microseconds to milliseconds after the big 
bang. Here, the only method to obtain information about the quark matter directly from the first 
principles of QCD is to numerically study QCD based on Monte Carlo simulations on the lattice. In 
the study of lattice QCD with dynamical quarks, the Boltzmann weight is proportional to the quark 
determinant. At nonzero chemical potentials, however, the quark determinant becomes complex 
and causes a serious problem - the Monte Carlo procedure is not justified because of the complex 
Boltzmann weight. In the low density region, because the fluctuations of the complex phase are 
small, we can avoid the problem by the reweighting method in which the complex phase is treated 
as a correction factor of the observables (reweighting factor). Using the reweighting method, we 
can also vary coupling parameters of the system by absorbing the difference of the Boltzmann 
weight at different coupling parameters into the reweighting factor. This is powerful in a study of 
the phase structure in which a survey over a range of coupling parameter space is mandatory. 

At higher densities, however, larger fluctuations of the complex phase introduce two severe 
difficulties. One is the sign problem. Because of large fluctuations of the complex phase in the 
reweighting factor, exponentially large statistics is required to obtain reliable estimates for the 
reweighted observables. Several methods have been proposed to remedy or mitigate the sign prob¬ 
lem py, 2). Another problem is the overlap problem. When we try to shift simulation parameters 
largely by the reweighting method, the reweighting factor tries to enhance part of the Boltzmann 
weight whose statistical quality is low. This can be easily seen by viewing the Boltzmann factor 
in terms of the histogram for relevant observables. When expectation values of observables vary 
largely with the shift of the simulation parameters, the reweighting factor has to enhance part of 
the histogram far from the original peak position. Because it is statistically quite hard to achieve 
highly accurate details of the histogram around such a point, the reweighting method may lead 
to completely unreliable results (see, e.g. ref. [3]). This makes it difficult to study high density 
QCD, in which the transition is expected to be of first order and thus expectation values can jump 
largely around the transition point. 

In this paper, we focus on the overlap problem. The overlap problem is expected to be milder 
if one changes a couple of parameters at the same time. In an early trial to identify the critical 
point in a high density region of QCD jJS], Fodor and Katz shifted the chemical potential and 
the gauge coupling (temperature) simultaneously along the crossover curve to achieve a better 
overlap (multiparameter reweighting method). More recently, the WHOT-QCD Collaboration 
investigated the phase structure of Af-flavor QCD in the heavy-quark region and found that the 
system at large quark masses is controlled by only two combinations of parameters, /3 + 48 ]C/^=i K / 
and KjI f t cosh(/ij/T), where (3 = 6 / g 2 is the gauge coupling, kj and [if are the hopping 

parameter and chemical potential for the / th flavor, and A/ is the temporal lattice size [3]. This 
means that, when one changes the coupling parameters while keeping these combinations constant, 
the system does not change and thus the overlap problem does not arise. We expect that similar 
combinations of parameters also exist in the light-quark region. 

In the study of ref. [3|], the multipoint reweighting method (?) for f3 played an important role: 
Combining configurations obtained at different /3, we could calculate the effective potential in a 
wide range of the observable values, which was mandatory in a reliable evaluation of the transition 
point. 

In the present paper, we extend the multipoint reweighting method to the multiparameter space 
of /3 and Kf, and test if the method helps to overcome the overlap problem in the light-quark region, 
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performing simulations in a simpler case of zero-density QCD. We measure histograms of physical 
quantities at a series of simulation points in two-flavor QCD and apply the multipoint reweighting 
method to calculate the meson masses as continuous functions of (3 and n. We then determine 
lines of constant physics in the (/3, k ) space and evaluate the derivatives of the lattice spacing with 
respect to (3 and n along the lines of constant physics (inverse of the beta functions), which are 
needed in a calculation of the equation of state at finite temperature. 

In the next section, the multipoint reweighting method is introduced, and in Sec lIIIl we examine 
the overlap problem by performing numerical simulations in two-flavor QCD. We then calculate the 
meson masses, the lines of constant physics, and the derivatives of the lattice spacing with respect 
to f3 and k along the lines of constant physics in Sec. IIVI Section [V] contains our conclusions. 


II. MULTIPOINT REWEIGHTING METHOD 
A. Multiparameter reweighting method 

Let us consider QCD with Nf flavors of quarks and define a histogram for a set of physical 
quantities X = ( X \, X- 2 , • • ■) by 


r Ni 

w(X;/3,k,h) = VU Y[S(Xi - Xi) e~ §G JJ det M(/3, K f , n } ). (1) 

J i f =1 

where Sq is the gauge action, M is the kernel matrix of the quark action, and X = {X\, X 2 , • • •) are 
the operators for X. The coupling parameters of the theory are the gauge coupling (3, the hopping 
parameters (ki,K 2 ,"‘ i K N t ), and the chemical potentials (hi, fi 2 , ■ ■ ■ ,A 4 N t )■ For simplicity, we 
denote the set of coupling parameters (/3, n 1 , • • • , /xi, • • ■) as b. 

Then, the partition function is given by Z(b) = fw(X; b) dX with dX = W^dXi, and the 
probability distribution function of X is given by Z _1 (6) w(X;b). The expectation value of an 
operator 0[X], which is written in terms of the operators X is evaluated as 

(0[X]) {b) = ^- ) f 0[X\ w(X; b) dX. (2) 

For convenience, we also define the effective potential as 

V eS (X-,b) =-In w(X-b). (3) 

Let us consider a calculation of the histogram at b using configurations generated at bo- Such 
a calculation can be easily done with the reweighting method by choosing S(b) and S(bo) as the 
first two elements of X, where 


S = S g -J2 In det M (4) 

/ 

is the effective action of QCD and S(b ) is the value of the action with the coupling parameters 
b evaluated on the configuration generated at the simulation point. Let us denote S(b) = S and 
S(bo) = S 0 , and redefine X as the set of remaining elements of X other than S and So- The 
histogram obtained by the simulation at bo is given by w(X, S , So; bo). From (P), we find that the 
histogram at b is simply given by 

w(X, S , So; b ) = w(X, S, S 0 ; bo), (5) 
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and the histogram of X at b is given by 


w{X- b) = J w{X, S , S 0 ; b) dS dS 0 . (6) 

In a simple case of quenched QCD, S = —6IV s ite/LP with P the plaquette and lV s ite the number of 
sites (lattice volume). Then, because both S and So are fixed when P is fixed, we have w(P] j3o) = 
w(P,S,S 0 -,fa), and 

w(P~fa = e 6Nsit ^~ MP w{P; fa), 

V eS (P ; fa = V e s(P; fa) - 6 N site (f3 - fa)P. 

If we approximate w(P',fa) by a Gaussian distribution centered at Po = (P)(b 0 ), he., w(P; fa) = 
exp[— V e s(P] /?o)] oc exp[— a(P — Po) 2 ] with an appropriate constant a, then the expectation value 
of P at (5 is given by P = Po + 3A^ s i te (/3 — fa) /a. When /3 — fa is large, P leaves the statistically 
reliable region of the original histogram rc(P; fa) and thus the results at f3 become unreliable (the 
overlap problem). As mentioned in the Introduction, in order to study the expected first-order 
transition of QCD at high densities, we need to obtain w and V e s reliably in a wide range of X. 


B. Multipoint multiparameter reweighting method 

To overcome the overlap problem, we extend the reweighting formulas to combine configurations 
obtained at different simulation points B0 for the case with dynamical quarks. We perform a 
series of N sp simulations at bj with the number of configurations N t where i = 1, ■ ■ • , N sp . Let us 
denote S = (S i, • • • , S/v sp ) with Si = S(bi). Using (JSJ), the probability distribution function at 6* 
is related to the histogram at b as 

Z~\bi) w(X, S, 5; bi ) = Z-\bi) e-( Si - s ) w(X, S, 5; b) (7) 

where S = S(b). We then obtain 

N sp 7V sp 

2 fa Z~\bi) w(X, S, S; bi) =e s Y J fa z ~\bi) e~ Si w{X , S, 5; b). (8) 

i —1 i=l 

Note that the left-hand side of © gives the naive histogram using all the configurations disregarding 
the difference in the simulation parameters From this relation, we find 

N s p 

w(X, S , 5; b ) = G(S, S; 6, b) X fa z ~\ h i) w (fa S i fa h) (9) 

1=1 


where 6 = (6i, • - - , b 7 v sp ) and 


G(S, S; 6,6) 


5 


T.tlfa^z-fabi) 


( 10 ) 


This expression means that the histogram w(X, S, S; b) at b is given by multiplying G(S, S;b,b) by 
the naive histogram. 
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To calculate G(S, S',b,b), we need the values of Z(bi), the partition function at bi. Here, we 
note that the partition function Z at (3 is given by 


Z(b) = / w(X, S, 5; b ) dXdSdS 


Wep - 

YNi G(S, s-,b,b) Z~\bi) w(X, S, S]bi) dXdSdS 


%— 1 
iVsi 


Y N i{G(S,S-,b,b) 


i=l 


(bi) 


( 11 ) 


which is just the naive sum of G(S, S; b, b) over all the configurations disregarding the difference in 
the simulation parameters. Then, Z(bi) can be determined by the consistency relations, 


N s , 


N sx 


Z(b i ) = Y N k{G(S,S;b i ,b )} = £> fc 


,-Si 


^ ti n \Z%We- B 'Z-Hbj) 


( 12 ) 


k=1 'fc=i 


for i = 1, ■ ■ • , AT sp , but up to an overall factor. Denoting /,; = — lnZ(6j), these equations can be 
rewritten as 


1 = Y Nk { I Y ex P[^* ~ 4' ~ fi + f: 

k= 1 \ \j=l 


jj | ) • * — I)'"') A^sp- 

' (b k ) 


(13) 


Starting from appropriate initial values of /), we solve (11311 numerically by an iterative method. 
Note that one of the /)’s must be fixed to remove the ambiguity corresponding to the undetermined 
overall factor. 

The expectation value of an operator X at b can be evaluated as 


N m 


(x ) {b ) = Xw(X,S,S-b)dSdS = Y(^Y N i{XG(S,S-,b,b))^ )} (14) 


which is just the naive sum of XG over all the configurations disregarding the difference in the 
simulation parameters. From this formula, we see that the histogram of X at b is given by 




w(X-,b) = Y N *( 6 (X-X)G(sJ-b,b)\ 

' ' (pi) 


1=1 


(15) 


III. TEST STUDY OF OVERLAP PROBLEM 
A. Two-flavor QCD 

To test the multipoint reweighting method, we perform simulations of QCD with degenerate 
two-flavor clover-improved Wilson quarks and RG-improved Iwasaki glues at zero density. The 
gauge action is given by 


S G = —6AQ t e /3 P, 


(16) 
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FIG. 1: The average of Slndet M/dn and its cubic spline interpolations on each configuration generated at 
P = 1.80 and k = 0.1400 (circle), 0.1425 (square) and 0.1440 (triangle). Statistical errors are estimated to 
be around the thickness of the curves. 


where N s a e = x A} is the lattice volume, and P is the improved plaquette by Iwasaki, 

P = c 0 W {lxl) + 2 Cl W^ lx2) , (17) 

with ci = —0.331, Co = 1 — 8ci [§j], and W ix i is the (i x j ) Wilson loop. The quark action is given 
by 

2 

Sq = ( 18 ) 

/=i x ,y 

Mxy = &xy ^ ^ ^ ^ (1 Oju) Px,fi^x+p,,y T (1 T 'J/i) ^xy ^SW ^ ^ ^ &[UjPxfiv 

\i y.>u 

where k is the hopping parameter common to two flavors, and F Xflu is the standard clover-shaped 
lattice field strength. For the clover coefficient cswi we adopt a mean-field value by substituting 
the one-loop result for the plaquette, csw = (1 — O.Sdlll/T -1 ) -3 / 4 . The improvement parameters of 
the action are the same as those adopted in refs. a S-Q. 

The simulations are carried out on an 8 4 lattice at 9 simulation points (all the combinations of 
P = {1.800, 1.825, 1.850} and k = {0.1400, 0.1425, 0.1440}) for the test study in this section, and 
on a 16 4 lattice at 30 points (all the combinations of P = {1.806, 1.8125, 1.819, 1.825, 1.831, 1.837} 
and k = {0.14000, 0.14125, 0.14250, 0.14300, 0.14400}) for the determination of lines of constant 
physics and beta functions in Sec. IIV1 The number of configurations for the measurement is 200 
at each simulation point, and statistical errors are estimated by a jackknife method. 


B. Reweighting with dynamical quarks 

The effective action S = Sq — l n det M consists of the gauge part Sq and the quark part 
IndetM. Calculation of the latter requires large computational cost. In this study, we evaluate it 
by measuring the first and second k derivatives of lndet M at several Ki s, where i = 1, 2, • • •, on 
each configuration, and interpolate ln det M between Ki and k.;+i assuming a quadratic function in 
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FIG. 2: Effect of the /9-dependence in csw on the value of 1 x 1 Wilson loop at k = 0.140. Two filled 
triangles represent the results obtained directly at the simulation points /? = 1.825 and 1.850. The green 
and red curves, which are almost completely overlapping with each other, are the results of the reweighting 
with and without the linear term in (1231) . 
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terms of k 4 , 

lndet M(k) = lndet M(«j) + C\{k — Ki) + — Ki) 2 + Cs(k — Ki) 3 + C±(k — Ki) 4 . (19) 


Then, the derivatives are written as 


d In det M 
8k 

8 2 In det M 

8 .^ 


(«) = Ci + 2 C 2 (k - Ki) + 3C 3 (k - Ki) 2 + 4C 4 (/s - At*) 3 

(k) = 2C*2 + 6C3(k — Ki) + 12C < 4(k — Ki) 2 . 


( 20 ) 

( 21 ) 


We fix the four coefficients C a such that the first and second derivatives reproduce the measured 
values at Ki and Ki + 1 , i.e., 


C, = 4b c 2 = I4 2) , c 3 = 


4 « + 2 4 2) 


44 , - 4 1 ’ 


3h 


d (2) +d (2) 

, C 4 = - ^ +1 ntq * + dt+1 A ^ 1 ■ (22) 


2 h 3 


4 h 2 


where d ^ = [<9In det Af/<9re]( k,), d^ 2> = [d 2 lndet M/dK 2 }(Ki) and h = «j + i — «i. Using the 
coefficients C a , we obtain lndet M(ac) parameterized by (1191) in the range of Ki < k < Ki+\ . This 
determines In det M up to an overall constant that is redundant in the reweighting calculations. 

The derivatives of In det M are given by traces of some combination of M -1 and derivatives of 
M (see ref. [2] for explicit expressions). We compute these traces by the random noise method. To 
reduce errors due to finite number of noise vectors, N no i se , we use the random noise method only for 
the trace over spatial indices, and we calculate the traces over color and spinor indices exactly. As 
discussed in ref. [2], this procedure helps us reduce iV no j se , in particular, with Wilson-type quarks. 
In this study, we adopt IY no j se = 5. 


1 We note that In det M is an even function of k. However, the first derivative exists at nonvanishing k, and the 
expansion is possible assuming cancellation of odd-power terms in k by the higher order terms. Although we 
can consider an alternative expansion in terms of, e.g., (k 2 — a: 2 ), the difference is absorbed by the higher order 
terms, and the quality of the fit is not improved in the present case. Because the derivatives in k are directly 
related to measurable observables, we prefer the expansion m- 
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FIG. 3: Left: The expectation value of the improved plaquette P = c$W lxl + 2ciW lx2 at /3 = 1.825. Black 
dots are the expectation values obtained by the simulations at k = 0.1400, 0.1425, and 0.1440. Blue curve 
is the result of single-point reweighting using the configurations at ft = 0.140 only. Right: Red, green, blue, 
magenta and light blue curves are the probability distribution functions of P at various ft’s obtained by the 
single-point reweighting method using the configurations at ft = 0.140 only. Black dashed curves are the 
original probability distribution functions at the simulation points ft = 0.1425 and 0.1440. 


In Fig. |TJ we show an example of the interpolation for <9 In det M/ds. Open symbols are the 
averages of this derivative measured at /3 = 1.80 and ft = 0.1400 (circle), 0.1425 (square) and 
0.1440 (triangle). The curves (with error bars) are the results of interpolation using the data of 
<91ndetM/<9ft and d 2 In det M/dn 2 with k = 0.1400, 0.1425 and 0.1440 on each configuration. 

The reweighting in the P direction requires care because the quark kernel M is dependent 
on P through the clover coefficient csw hi our choice. We take into account the effect of the P 
dependence in csw by a linear approximation 


m 


In det M ( P , ft) 


In det M (/3 0 , ft) + (P — /3 0 ) 


dcsw d In det M 


dp 


dc S w 


Ad,« 


(23) 


In Fig. [2l we show the results of 1 x 1 Wilson loop at ft = 0.140 with and without the linear term 
in (|23l) . We find that the differences are at most 0.05% and are much smaller than the statistical 
errors. We thus consider that the effects of the P dependence in csw is small and thus the linear 
approximation is sufficiently safe in the range of the coupling parameters we study. In the following, 
we adopt the linear approximation ([23]) . 


C. Overlap problem and multipoint reweighting 

In the left panel of Fig. [3l we show the results for the improved plaquette P = coW lxl +2c\W lx2 
of the Iwasaki action at /3 = 1.825. The black dots represent the expectation values of P at the 
three simulation points without reweighting. The blue curve shows the results of the single-point 
reweighting method using the configurations at ft = 0.1400 only. We note that the blue curve fails 
to reproduce the data at ft = 0.1425 and 0.1440. Even the error bars based on a standard jackknife 
analysis are unreliable. The reason can be easily understood by consulting the histogram of P: The 
red curve in the right panel of Fig. [3] is the original probability distribution function at (/3, ft) = 
(1.825,0.140), and green, blue, magenta and light blue curves are the probability distribution 
functions at ft = 0.1412, 0.1425, 0.1434, and 0.1440, respectively, predicted by the single-point 





































FIG. 4: Left: The expectation value of P = cqW 1x1 + 2ciW lx2 as a function of k at /3 = 1.825. Black dots 
are the expectation values obtained by the simulations at ft = 0.1400, 0.1425, and 0.1440. Blue, green and 
purple curves are the results of the single-point reweighting method using the configurations at ft = 0.1400, 
0.1425, and 0.1440, respectively. Red curve is the result of the multipoint reweighting method combining the 
configurations at the three k’s. Right: Red, green, blue, magenta and light blue curves are the probability 
distribution functions of P by the multipoint reweighting method at various k’s. Black dashed curves are 
the original probability distribution functions at the three simulation points, ft = 0.1400, 0.1425, and 0.1440. 
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FIG. 5: The /3-dependence (left) and ft-dependence (right) of the histogram for N si ^ e (dS/dP) and 
N^[dS/dn] sub = N^[dS/dK - (2883Vf ft 4 /c 0 ) (9S 1 /<9/3)], where N site = 8 4 . 


reweighting ([5]) using the configurations at k = 0.1400. Because a probability distribution function 
at ft other than the simulation point is calculated as a product of the reweighting factor and the 
original distribution function, the shifted distribution functions cannot go out of the range of the 
original distribution P ~ 1.64-1.695 at ft = 0.1400, and they fail to reproduce the true distribution 
functions at ft = 0.1425 and 0.1440 shown by the black dot-dashed curves. Accordingly, the 
expectation value, which is approximately the peak position of the distribution function, cannot 
go out of the range of the original distribution, as shown in Fig. |4] (left). Furthermore, due to 
the poor statistics around the boundary of the original distribution, we cannot estimate the errors 
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FIG. 6: The histogram of N si ^ e (dS/df3) and N^JOS/dn) , normalized by the maximum height. 


there reliably. 

The multipoint reweighting method introduced in the previous section can enlarge the range of 
reliable reweighting by combining different ranges of distribution obtained at different simulation 
points. The result of multipoint reweighting combining the configurations at k, = 0.1400, 0.1425, 
and 0.1440 is shown by a red curve in the left panel of Fig. [4j Results of the single-point reweighting 
method using the configurations at k = 0.1400, 0.1425, and 0.1440 separately are also shown by 
blue, green and purple curves, respectively. We find that, unlike the case of single-point reweighting, 
the red curve smoothly connects all the simulation results with small errors. In the right panel 
of Fig. [U probability distribution functions from the multipoint reweighting method are plotted 
for k = 0.1412, 0.1424, 0.1436, and 0.144. The probability distribution functions reproduce and 
smoothly interpolate the original distribution functions at different simulation points. 

The method is applicable to other observables too. In the calculation of the equation of state, 
we calculate the following combination of energy density(e) and pressure (p), 


e — 3 p 

2^4 


= Nf 


= Nf 


1 


NfN t 

dp 


da 


dS_ 

da 

1 


o 

dS 


dn 

Nj*N t dp l n da 


1 dS 
N^N t ~dK 


o J 


(24) 


where (• • • )o is the expectation value at finite temperature with the zero-temperature value sub¬ 
tracted. We then need the expectation values of the derivatives of the action S = Sc — J2f lndet M, 
i.e. dS/dp and dS/dn. 

From our experience in the heavy quark region (l3 |. we expect that the n dependences of dS/dp 
and dS/dn are strongly correlated with each other. We thus consider the combination 


dS_ 

dn 


J SUB 


dS_ 

dn 


288AW dS 

w 


co 


(25) 


as a component approximately perpendicular to dS/dp, by subtracting the leading order contri¬ 
bution in the hopping parameter expansion from dS/dn. In the right and left panels of Fig. [5l 
we show the results of the multipoint reweighting method for the P and n dependence of the two- 
dimensional histogram of dS/dp and [(9S'/cIk]sub- The histograms are obtained by combining the 
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FIG. 7: The pseudoscalar meson mass raps& (left) and vector meson mass my a (right) as functions of 1 /ft. 
Data points without a clear plateau are removed. 
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FIG. 8: The pseudoscalar and vector meson mass FIG. 9: The lines of constant physics at mps/my = 
ratio mps/my as functions of 1/ft. 0.70, 0.72, 0.74, and 0.76 in the (/?, ft) plane. 


configurations at 9 simulation points (3 /3’s x 3 ft’s) on the 8 4 lattice. We see that the histogram 
moves smoothly as /? and ft are varied, without hitting a boundary. We can thus compute the 
expectation values of dS/d/3 and [9S'/9ft]suB as continuous functions of /? and ft in this range of 
the coupling parameters, without the overlap problem. 

The applicable range of the multipoint reweighting method can be estimated easily from the 
sum of the histograms measured at each simulation point. As explained in Sec. HU the expectation 
values (fT4|) in the multipoint reweighting method are obtained by just the naive sum of the operators 
multiplied by the reweighting factor G over all the configurations disregarding the difference of f3 
and ft. We show, in Fig. [6J the contour plot of the naive sum of the histograms obtained at all 
9 simulation points in the (dS/d/3, dS/dn) plane. In the region painted by bright colors, many 
configurations are available, and thus a reliable calculation is possible. 


IV. LINES OF CONSTANT PHYSICS AND BETA FUNCTIONS 

Because the multipoint reweighting method enables us to compute observables as continuous 
functions of /3 and ft, it is useful to determine the lines of constant physics in the coupling parameter 
space as well as the beta functions, a(d(3/da) and a(dK/da). These quantities are needed in a 
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FIG. 10: The /3-dependence of mpsa (left) and my a (right) on the line of constant physics at mps/my = 
0.70. Results of n-th order polynomial fits are also shown. 



FIG. 11: The same as Fig. [TOlbut at mps/my = 0.74. 


calculation of the equation of state. To calculate them, we perform simulations at 30 simulation 
points (all the combinations of (3 = {1.806, 1.8125, 1.819, 1.825, 1.831, 1.837} and k = {0.14000, 
0.14125, 0.14250, 0.14300, 0.14400}) on a 16 4 lattice. We then combine the configurations at the 
30 simulation points by the multipoint reweighting method. We determine IndetM as a function 
of k using the interpolation method discussed in Sec. IIIIB1 For this calculation, we compute the 
derivatives of In det M by the random noise method with the number of random vectors lV no i se = 5 
at the five k points. 

In this study, we define the lines of constant physics by fixing the dimensionless ratio of pseu¬ 
doscalar and vector meson masses, mps/my = mpsa/mya, in the (/3, k) space, where a is the 
lattice spacing. Along a line of constant physics thus defined, a varies as we change (3 or k. The 
beta functions are defined as the derivatives of (3 and k by a along lines of constant physics. 

To determine the meson masses as functions of [3 and k, we need meson correlation functions 
Q(t) at various f3 and k. Because the computational cost for these correlation functions is relatively 
low, we compute them at more points of f3 and k than the simulation points using the multipoint 
reweighting method. For j3 , we choose 31 points, inserting five additional points between each 
two succeeding simulation points at (3 = 1.806, 1.8125, 1.819, 1.825, 1.831, and 1.837. For k, we 
choose 43 points at k = 0.1400, 0.1401, • • •, 0.1442. (At k > 0.1442, we cannot always find a stable 
plateau in the effective mass plot discussed below.) At each measurement point, we measure Q(t) 
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FIG. 12: The K-dependence of mpsa (left) and my a (right) on the line of constant physics at Tops /my = 0.70. 
Results of n-th order polynomial fits are also shown. 




FIG. 13: The same as Fig. [I2]but at mps/my = 0.74. 


on 200 configurations every 10 trajectories after thermalization. We also average over 8 different 
source points. 

We calculate mpsa and my a by a cosh fit of the meson correlation function in the range t/a = 5- 
8, when a plateau of effective mass is identified in the range. The results of mpsa and mya are 
shown in the left and right panels of Fig. [TJ respectively. The errors are estimated by the jackknife 
method. The mass ratio mps/my is shown in Fig. [HJ 

At each (3, we interpolate m-ps/my as a function of k to determine the lines of constant physics 
for mps/my = 0.70, 0.72, 0.74 and 0.76. The results are shown in Fig. [9j Because n is more 
sensitive to mps/my than f3, the values of k for a line of constant physics greatly vary with 
mps/m-v- 

Along a line of constant physics, mpsa and mya vary with j3 or k. The (3 and k, dependences 
of these masses at mps/my = 0.70 and 0.74 are plotted in Figs. flOl and fill and Figs. fl2l and [T3l 
respectively. The results at mps/my = 0.72 and 0.76 are similar. We note that, although the 
errors estimated at each point vary, the central values form quite smooth curves. This may be due 
to the correlation among different points by the reweighting procedure. 

In Figs. [Toll 131 results of n th order polynomial fits are shown by dotted lines for n = 1-4. We 
find that the masses are well fitted with n = 1 or 2 in the range of (3 and k that we study. 

From the fit functions, we can calculate the derivatives, (mpsa) db/d(mpsa) and 
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FIG. 14: Beta functions a(d(3/da) (left) and a(da/da) (right) determined in terms of the quadratic fits 
(n = 2) of my a. The dashed lines are the one-loop perturbative values of a(df3/da) in QCD with zero and 
two flavors of massless quarks. 




FIG. 15: The same as Fig. [HI but by the linear fits (n = 1). 


(mya) db/d(mya) with b = /3 and k, along the lines of constant physics. Because both mps 
and my are constant on a line of constant physics, we expect 


(m PS a) db 

a(mpso) 


db db 

{mva) d(^a) = 


b = ft, n. 


We confirm that the beta functions in terms of mps« and my a are almost indistinguishable from 
each other: The differences are at most 0.15% in the range of coupling parameters that we study 
and are at most 14% of the statistical errors. In the following, we adopt mya for the scale. 

Consulting Figs. fT0l[T3l we adopt the quadratic fits (n = 2) of mya for the calculation of the 
beta functions a(d/3/da) and a(dn/da ). The results are shown in Fig. fl4l The errors shown are 
statistical only. Recall that, because the data at different coupling parameters are correlated due to 
the reweighting procedure, the statistical errors of the beta functions turn out to be much smaller 
than those from a naive impression of the meson mass plots. To get an idea about the magnitude 
of systematic errors due to the fit ansatz of mya, we also show the results with linear (n = 1) 
ansatz in Fig. [T5l 

In the left panel of Fig. ITT1 the results of one-loop perturbation theory for a(dfd/da) with zero 
and two flavors of massless quarks are shown by dot-dashed and dashed lines, respectively. Taking 
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into account the systematic errors, our results are approximately consistent with the two-flavor 
perturbative value. On the other hand, we expect that a(dn/da ) approaches zero in the large /3 
limit. Such tendency is not visible yet in the right panel of Fig. Q3J 


V. CONCLUSIONS AND OUTLOOK 


We studied the multipoint reweighting method in a multidimensional parameter space to avoid 
the overlap problem. Performing simulations in two-flavor QCD with Iwasaki’s improved gauge 
action and improved clover quark action, we find that the overlap problem can be avoided by 
appropriately combining configurations at different simulation points by the multipoint reweighting 
method. We have further shown that the method is useful in calculating the line of constant physics 
as well as the beta functions, which are required in the evaluation of thermodynamic properties such 
as the equation of state. Extending the multipoint reweighting method to the reweighting study 
on anisotropic lattices 14|, we can also calculate the Karsch coefficients [15j which are required in 
the evaluation of the equation of state by the differential method . 

Our final objective is to carry out a study of finite density QCD using the multipoint reweighting 
method. In our previous study in the heavy quark region [3,1111], we found that the leading effects of 
the chemical potential can be absorbed by a shift of coupling parameters and the overlap problem 
is avoided by keeping these shifted parameters constant. We expect that a similar shift to absorb 
the main effects of the chemical potential also exists at lighter quark masses. Combining with the 
multipoint reweighting method, we may be able to investigate the phase structure of finite density 
QCD at lighter quark masses, maximally avoiding the sign problem. 
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